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Subgrid scale modeling in solar convection 
simulations using the ASH code 

By Y.-N. Young, M. Mieschf and N. N. Mansour 


1. Motivation and objectives 

The turbulent solar convection zone has remained one of the most challenging and 
important subjects in physics. Understanding the complex dynamics in the solar con- 
vection zone is crucial for gaining insight into the solar dynamo problem. Many solar 
observatories have generated revealing data with great details of large scale motions in 
the solar convection zone. For example, a strong differential rotation is observed: the 
angular rotation is observed to be faster at the equator than near the poles not only 
near the solar surface, but also deep in the convection zone. On the other hand, due 
to the wide range of dynamical scales of turbulence in the solar convection zone, both 
theory and simulation have limited success. Thus, cutting edge solar models and numer- 
ical simulations of the solar convection zone have focused more narrowly on a few key 
features of the solar convection zone, such as the time-averaged differential rotation. For 
example, Brun & Toomre (2002) report computational finding of differential rotation in 
an anelastic model for solar convection. A critical shortcoming in this model is that the 
viscous dissipation is based on application of mixing length theory to stellar dynamics 
with some ad hoc parameter tuning. 

The goal of our work is to implement the subgrid scale model developed at CTR into 
the solar simulation code and examine how the differential rotation will be affected as a 
result. Specifically, we implement a Smagorinsky-Lilly subgrid scale model into the ASH 
(anelastic spherical harmonic) code developed over the years by various authors. Readers 
are referred to (Clune et al. 1999) and (Miesch 1998) for a detailed description of the 
ASH code. This paper is organized as follows. In §2 we briefly formulate the anelastic 
system that describes the solar convection. In §3 we formulate the Smagorinsky-Lilly 
subgrid scale model for unstably stratified convection. We then present some preliminary 
results in §4, where we also provide some conclusions and future directions. 


2. The anelastic equations of motions 

In the solar convection zone, the depth extends over many density scale heights and 
thus the effects of stratification need to be captured despite the fact that acoustic 
timescales are much shorter than the large-scale convection timescales. An appropri- 
ate approach is to filter out sound waves while maintaining the density stratification, 
namely the anelastic approximation first proposed by Gough (1969) and later adapted to 
the solar convection zone by Gilman & Glatzmaier (1981). In stellar models, a fluid layer 
is unstable if the layer is superadiabatic. The degree of superadibaticity is quantified as 
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(2.1) 




C P 


where d is the depth of the convection zone, Cp the specific heat capacity at constant 
pressure, s the entropy, and r the radial coordinate. The entropy gradient is calculated 
at some fiducial level in the convection zone. In the sun, the departure from adiabaticity 
in the convection zone is extremely small (Christensen-Dalsgaarcl et al. 1993), 


e < 10 


-4 


( 2 . 2 ) 


If we assume that the velocity in the horizontal (transverse to gravity) directions is mostly 
driven by pressure gradients, the coupling between vertical and horizontal velocities leads 
to the following relationship between the convection velocity v and sound speed c s 

1, (2.3) 

c s 

where A4 is the Mach number of the flow. 

The basic idea in the anelastic approximation is to expand the variables in terms of 
the small parameter e, and collect terms of the zeroth and first orders. As is usually 
found in perturbation theory, the zeroth order equations describe the basic state, which 
remain stationary over timescales at which the first order variables vary. The fundamental 
equations are the compressible, stratified Navier-Stokes equations in a rotating reference 
frame 


p + (v • V)v^j = -VP + pg + 2v x fl 0 + V • V, 
ds 

p&tt + P©v • Vs = —V 

ot 

S + V ■ (pv) = 0, 


1 q + 2 vp (eijeij - i(V • v) 2 ^j - pE, 


(2.4) 

(2.5) 

( 2 . 6 ) 


where s is the entropy, 0 is the temperature, and E is the nuclear energy generation rate. 
The viscous stress T> is defined as 


V = 2 pv - * (V • v)6ij i , 

and q is the composite non-convective heat flux, defined as 

q = —KpQ\7s — K r pCpS/Q. 


(2.7) 


( 2 . 8 ) 


and the various transport coefficient such as the viscosity v, thermal diffusion k, and 
radiative thermal diffusion n r are assumed to be functions of the radial coordinate only. 

Upon substituting the variables expanded in terms of the small parameter e, at ze- 
roth order we obtain the hydrostatic equation. Denoting the zeroth order terms with a 
subscript 0, the equations at the leading order are 


dPo 

dr 


+ A = -p 0 g, 


Po = PpoBo, 


(2.9) 

(2.10) 
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A = ([p 0 (v • V)v - V • V - 2p 0 (v x f2 0 )] r ) t ■ (2.11) 

The angular brackets in Equation (2.11) denote horizontal averages, and the square 
brackets denote the radial component of the enclosed vector. We evaluate A only as we 
update the reference state (denoted with a subscript 0). 

At the next order, we obtain the anelastic equations for compressible, stratified fluids 


^ + p 0 (v • V)v = -VP + V • V + pg + 2p 0 (v x fi 0 ) + A f, (2.12) 
ds 

PoQ °~di + Po ®° v ' V ( s ° + s ) = v ' («Po@oV(s 0 + s) + K r p o C'pV(0 o + 0)) 

(V • v) 2 ) - Po £ , (2.13) 

V • (p 0 v) = 0, (2.14) 



and 


p _ P 0 _ P s 

Po Po 0o 7 Po Cp 


(2.15) 


We note that both the reference state and perturbation terms have been retained in 
Equation (2.13) because the gradient of the reference state entropy (so) varies in ampli- 
tude from O(e) in the convection zone to 0(1) in the convectively stable region. This 
scaling is also true for the term involving thermal diffusion of the reference temperature 
( K r d&o/dr ). We retain the term n r dQ/dr because we wish to allow for small deviations 
of the radiative heat flux from spherical symmetry. 

The viscosity coefficient v is assumed to be a function of the reference density po. The 
functional dependence is obtained from mixing length theory, and the specific value is 
chosen in an ad hoc fashion for a given numerical resolution and combination of parame- 
ters. Our main goal in this work is to replace these “mixing length” viscosity and thermal 
diffusivity with those based on Smagorinsky subgrid scale model. More details will be 
presented in §3. 

The reference state is crucial in the simulation of the anelastic system. In the ASH code, 
the reference state is updated only every few hundred steps of advancing the anelastic 
equations. During the update process, the following equations are solved simultaneously 
for the reference pressure Pq and density pq\ 


1 ds 0 
Cp~dr 

dPo 

dr 


1 d In P 0 
7 dr 

-Po9- 


d In p 0 
dr 


(2.16) 

(2.17) 


Using the reference entropy gradient ( dso/dr ) and gravity (g(r)) profiles from the 
previous time step, Equations (2.16) and (2.17) are solved for the reference pressure 
and density. Figure 1 is an example from one of the ASH simulations. With similar 
reference states and viscosity profiles to those illustrated in figure 1, simulations of the 
ASH code have successfully generated differential rotation profiles that are similar to 
the solar observations. An example from ASH simulations is illustrated in figure 2. The 
numerical resolution for figure 2 is 65 x 128 x 256 in the radial, latitudinal, and azimuthal 
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Figure 1 . Reference state from ASH simulations, the horizontal axis is radial coordinate in units 
of solar radius. Panel (a): Reference pressure (dyne/cm 2 ) within the convection zone. Panel (b): 
Reference density (g/cm 3 ). Panel (c): the “mixing length” viscosity (cm 2 /sec) as a function of 
r. 


directions. The left panel is a filled contour plot of the averaged angular velocity as a 
function of latitude and radius, while the right panel shows the angular velocity at five 
latitudes as a function of radius from the bottom to the top of the convection zone. For 
this particular case, the viscosity is a prescribed function of radius as shown in figure 
1(c), and the Prandtl number is fixed at 1/4. This corresponds to the AB case in Brun 
& Toomre (2002), in which the differential rotation inside the convection zone is most 
prominent among all the simulation data presented in Brun & Toomre (2002). Figure 
2 is from simulations conducted on the NAS A/ Ames SGI cluster with different initial 
conditions. A typical differential rotation from solar helioseismic observations is shown in 
figure 2(c) (same as figure 1(b) in Brun & Toomre (2002)). The close similarity between 
observational solar differential rotation and the simulation data shows great promise that 
better results may be obtainable if the conveciton model is refined by building in more 
physics. For example, the viscosity and thermal diffusivity may be replaced by better 
dynamical models. Addition of a mechanism to couple the convection zone with the 
tachocline below may result in better agreement between simulation and observational 
data. As a preliminary initiative we have undertaken is to first improve the model for 
viscosity using Smagorinsky-Lilly’s model. 
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Figure 2. Panels (a) and (b): Differential rotation for the AB case in Bran & Toomre (2002) 
with different initial conditions. Panel (c) is the time-averaged rotation rates from five years of 
GONG data. 


3. Smagorinsky-Lilly subgrid scale model 

The principle underlying the Smagorinsky subgrid scale model is the balance of pro- 
duction of subgrid-scale turbulent kinetic energy and dissipation of isotropic turbulence 
energy at the characteristic eddy size. Let the subgrid turbulent production rate be de- 
noted by V, then for a given turbulent viscosity vt, the subgrid turbulent production 
rate is 

V = 2v T e ij e ij . (3.1) 

The dissipation at scale l is ~ q 3 /l, where q is a characteristic turbulent velocity. Finally, 
to find an expression for vt in terms of the resolved quantities, we utilize the Prandtl’s 
assumption: vt = Ciql, and upon substituting q into the dissipation rate and equating 
dissipation and production rates, we obtain the Smagorinsky subgrid scale model 

v t = (i C s l) 2 yj2eijeij , (3.2) 


where C s is a constant in the range of 0.1 < C s < 0.3. Lilly (1962) later extended the 
Smagorinsky model to turbulence in unstably stratified convection. Due to the stratifi- 
cation an additional parameter, the Richardson number, appears in the expression for 
the eddy viscosity. The idea that led to Lilly’s eddy viscosity is very similar to the above 
argument underlying the Smagorinsky model: buoyancy production or consumption must 
also be accounted for in the subgrid energy balance. Thus in the stratified case, the dis- 
sipation consists of contribution from both eddy dissipation and thermal diffusion. The 
Smagorinsky-Lilly eddy viscosity is simply the modified Smagorinsky model with the 
inclusion of stratification effect 


vt 


(C s Z) 2 \/2eyej 


i - 

vt 


n 0.5 


(3.3) 


Here kt is the eddy thermal diffusivity, and the flux Richardson number is defined as 


. _ g d In 6 
2e, 5 ejj dr 


(3.4) 
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Figure 3. The left pannel is the reference entropy gradient (erg/g K m) and the right pannel 
is the averaged Smagorinsky- Lilly viscosity (cm 2 /sec) as functions of r in units of solar radius. 


Alternatively, the flux Richardson number can also be defined through entropy as 

. _ g dsg/dr 


Rif = 


2c v y tj'j 


where dso/dr is the reference entropy gradient. In spherical geometry, we choose the 
length scale l as 


where A(r) is the radial grid spacing at location r, L max is the maximum (latitudinal) 
angular mode index. 

Figure 3 shows the reference entropy gradient and the averaged Smagorinsky-Lilly 
viscosity as a function of radius. Near the bottom of the convection zone, the reference 
entropy gradient is positive, implying a stably stratified region. This corresponds to the 
small bump seen in the Smagorinsky-Lilly viscosity profile. 

In the simulations we start from a stationary configuration with a negative reference 
entropy gradient throughout the convection zone. The boundary conditions for the ve- 
locity are stress free, which are not necessarily realistic but numerically convenient. As 
we are mostly interested in the differential rotation profile in the convection zone, the 
parameters and initial reference state are chosen to be the same as the AB case in Brun 
& Toomre (2002). We first obtain a turbulent convection zone using the mixing length 
viscosity and thermal diffusivity. As the simulation progresses, the differential rotation 
profile reaches a statistically stationary state. We then replace the mixing length vis- 
cosity and thermal diffusivity with the Smagorinsky-Lilly model with some prescribed 
constant coefficient C s . At present we average Smagorinsky-Lilly’s eddy viscosity over 
latitudinal and azimuthal angles. Strictly speaking this is inconsistent with the energy 
balance principle that gives rise to the eddy viscosity. We are now implementing a local, 
three-dimensional version of the Smagorinsky-Lilly model, which requires a considerably 
tedious alteration of the ASH code. For the following results, we present ASH simula- 
tions using the averaged Smagorinsky-Lilly eddy viscosity as a test study to see how the 
differential rotation may differ from that obtained with the “mixing length” viscosity 
and thermal diffusivity. 


A(r)r 2 

I J max(I J max + i) 
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Figure 4. Eddy viscosity (cm 2 /sec) using the length scale defined in Equation (4.1) (solid line) 
and the mixing length viscosity (dashed line). On the right is the functional form of f(r) used 
in constructing the eddy viscosity. 


4. Preliminary results and future directions 

Our experiences so far show that the averaged Smagorinsky-Lilly eddy viscosity (as 
shown in figure 3) is too small, and simulations (with resolutions 65 x 128 x 256) blow up 
soon after the Smagorinsky-Lilly viscosity is activated. A significant difference between 
the mixing length viscosity and the averaged Smagorinsky-Lilly viscosity is that the eddy 
viscosity is much smaller at the boundaries. Similar issues arise in large eddy simulations 
of the atmospheric flows, and it is usually necessary to incorporate a smooth transition 
from the eddy viscosity to a wall-model viscosity near the boundary. In the solar case, we 
can combine the mixing length viscosity with the eddy viscosity using a choice of length 
scale l as follows 

l 2 = fWsL + (1 - (4.1) 

where Isl is the length scale in Equation (3.6) and Im is defined as 


0.5 


Im = w 


Fm 


Cs l sj2eijeijyjl - Rif 


(4.2) 


Here /(r) consists of two error functions so that /(r) = 1 in the bulk of the convection 
zone, and /(r) = 0 at the boundaries of the convection zone. Figure 4 shows the resul- 
tant eddy viscosity (left panel) for a given choice for /(r) (right panel). Figure 5 shows 
the differential rotation at a instant in time from simulations using the length scale in 
Equation (4.2). 

It is clear from comparing figure 5 and figure 2 that the differential rotation rates 
are very different near the top of convection zone. In fact the solar observation suggests 
that the rotation rates decrease near the top of convection zone, which is the opposite 
of the trend in figure 5. However, near the bottom and in the bulk of the convection 
zone, the rotation rates from the ASH simulations using the combined eddy viscosity 
are very similar to those in figure 2. This implies that we may expect to find similar 
differential rotation profile once we implement a better model than Equation (4.2) for 
the eddy viscosity near the top of convection zone. Better results may be obtained as 
well once we implement the fully local Smagorinsky-Lilly’s eddy viscosity. 
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Figure 5. Panels (a) and (b): Differential rotation using the eddy viscosity with a transition 
to the mixing length viscosity near the boundaries of convection zone. 
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